function [ResEuler,yGrid] = EulerEqErrorPruning_mp(model,orderApp,setup)
 
% Some indices 
ny      = model.ny;        %Number of output variables
nx      = model.nx;        %Number of state variables 
derivs  = model.derivs;    %Deviatives as stored by Levintal 
 
%Unfold params
BETTA= model.params.BETTA;
B= model.params.B;
CHI= model.params.CHI;
CHI0= model.params.CHI0;
THETA= model.params.THETA;
DELTA= model.params.DELTA;
ALFA= model.params.ALFA;
PHI= model.params.PHI;
PHIzero= model.params.PHIzero;
ZI= model.params.ZI;
KAPAw= model.params.KAPAw;
ETA= model.params.ETA;
RHOR= model.params.RHOR;
PHIpai= model.params.PHIpai;
PHIpai_1= model.params.PHIpai_1;
PHIy= model.params.PHIy;
PHIy_1= model.params.PHIy_1;
PHIc= model.params.PHIc;
PHIc_1= model.params.PHIc_1;
PHIl= model.params.PHIl;
NU= model.params.NU;
U0= model.params.U0;
U0d= model.params.U0d;
OMEGAz= model.params.OMEGAz;
OMEGAd= model.params.OMEGAd;
OMEGAn= model.params.OMEGAn;
OMEGAp= model.params.OMEGAp;
OMEGAa= model.params.OMEGAa;
RHOz= model.params.RHOz;
RHOd= model.params.RHOd;
RHOn= model.params.RHOn;
RHOp= model.params.RHOp;
RHOa= model.params.RHOa;
Kss= model.params.Kss;
OUTPUTss= model.params.OUTPUTss;
PAIss= model.params.PAIss;
MUZss= model.params.MUZss;
Css= model.params.Css;
AA= model.params.AA;
lss= model.params.lss;
Rss= model.params.Rss;
Wss= model.params.Wss;
 
% The raw grid for the states 
xGrid_f = setup.xSim_decom.xt_f; 
if orderApp >= 2 
    xGrid_s = setup.xSim_decom.xt_s; 
else 
    xGrid_s = setup.xSim_decom.xt_f*NaN; 
end 
if orderApp >= 3 
    xGrid_rd = setup.xSim_decom.xt_rd; 
else 
    xGrid_rd = setup.xSim_decom.xt_f*NaN; 
end 
if orderApp >= 4 
    xGrid_4th = setup.xSim_decom.xt_4th; 
else 
    xGrid_4th = setup.xSim_decom.xt_f*NaN; 
end
if orderApp >= 5 
    xGrid_5th = setup.xSim_decom.xt_5th; 
else 
    xGrid_5th = setup.xSim_decom.xt_f*NaN; 
end
 
 
% The grid for numerical itegration 
tmp_e        = cell(1,model.ne); 
tmp_w        = cell(1,model.ne); 
for i=1:model.ne 
    tmp_e(i) = mat2cell(setup.GH_e,length(setup.GH_e),1); 
    tmp_w(i) = mat2cell(setup.GH_w,length(setup.GH_w),1); 
end 
eGrid = cartprod(tmp_e); 
wGrid = prod(cartprod(tmp_w),2); 
 
%% Building the residuals 
ResEuler = zeros(size(xGrid_f,2),ny+nx); 
yGrid    = zeros(size(xGrid_f,2),ny); 
%start of multiprocessing
numCPUs = feature('numcores');
parpool('local',numCPUs);
parfor (k=1:size(xGrid_f,2),numCPUs)
    % Setting the dimensions for the state components 
    xf_cu    = [zeros(nx,1);1]; 
    xs_cu    = [zeros(nx,1);0]; 
    xrd_cu   = [zeros(nx,1);0]; 
    x4th_cu  = [zeros(nx,1);0]; 
    x5th_cu  = [zeros(nx,1);0]; 
 
    xf_cup   = [zeros(nx,1);1]; 
    xs_cup   = [zeros(nx,1);0]; 
    xrd_cup  = [zeros(nx,1);0]; 
    x4th_cup = [zeros(nx,1);0]; 
    x5th_cup = [zeros(nx,1);0]; 
 
    xf_cu    = xGrid_f(:,k); 
    if orderApp>=2 
        xs_cu  = xGrid_s(:,k); 
    end 
    if orderApp>=3 
        xrd_cu  = xGrid_rd(:,k); 
    end 
    if orderApp>=4 
        x4th_cu  = xGrid_4th(:,k); 
    end 
    if orderApp>=5 
        x5th_cu  = xGrid_5th(:,k); 
    end 
 
    % Getting the states at time t 
    c_ba1 = model.h0(1,1) + xf_cu(1,1) + xs_cu(1,1) + xrd_cu(1,1) + x4th_cu(1,1) + x5th_cu(1,1);
    muz_cu = model.h0(2,1) + xf_cu(2,1) + xs_cu(2,1) + xrd_cu(2,1) + x4th_cu(2,1) + x5th_cu(2,1);
    d_cu = model.h0(3,1) + xf_cu(3,1) + xs_cu(3,1) + xrd_cu(3,1) + x4th_cu(3,1) + x5th_cu(3,1);
    n_cu = model.h0(4,1) + xf_cu(4,1) + xs_cu(4,1) + xrd_cu(4,1) + x4th_cu(4,1) + x5th_cu(4,1);
    paistar_cu = model.h0(5,1) + xf_cu(5,1) + xs_cu(5,1) + xrd_cu(5,1) + x4th_cu(5,1) + x5th_cu(5,1);
    a_cu = model.h0(6,1) + xf_cu(6,1) + xs_cu(6,1) + xrd_cu(6,1) + x4th_cu(6,1) + x5th_cu(6,1);
 
    % Getting the states at time t+1 
    xf_cup(1:end-1,1)     = derivs.hx*xf_cu; 
    if orderApp>=2 
        xf2_cu            = kron(xf_cu,xf_cu); 
        xs_cup(1:end-1,1) = derivs.hx*xs_cu + derivs.hxx*xf2_cu/2; 
    end 
    if orderApp>=3 
        xf3_cu             = kron(xf2_cu,xf_cu); 
        xf_xs_cu           = kron(xf_cu,xs_cu); 
        xrd_cup(1:end-1,1) = derivs.hx*xrd_cu + derivs.hxx*(2*xf_xs_cu)/2 ... 
                            + derivs.hxxx*xf3_cu/6; 
    end 
    if orderApp>=4 
        xf4_cu             = kron(xf3_cu,xf_cu); 
        xf2_xs_cu          = kron(xf_cu,xf_xs_cu); 
        xs2_cu             = kron(xs_cu,xs_cu); 
        xf_xrd_cu          = kron(xf_cu,xrd_cu); 
        x4th_cup(1:end-1,1)= derivs.hx*x4th_cu + derivs.hxx*(2*xf_xrd_cu + xs2_cu)/2 ... 
                            + derivs.hxxx*(3*xf2_xs_cu)/6 ... 
                            + derivs.hxxxx*xf4_cu/24; 
    end 
    if orderApp>=5 
        xf5_cu             = kron(xf4_cu,xf_cu); 
        xf3_xs_cu          = kron(xf_cu,xf2_xs_cu); 
        xf_xs2_cu          = kron(xf_cu,xs2_cu); 
        xf2_xrd_cu         = kron(xf_cu,xf_xrd_cu); 
        xs_xrd_cu          = kron(xs_cu,xrd_cu); 
        xf_x4th_cu         = kron(xf_cu,x4th_cu); 
        x5th_cup(1:end-1,1)= derivs.hx*x5th_cu ... 
                            + derivs.hxx*(2*xf_x4th_cu+2*xs_xrd_cu)/2 ... 
                            + derivs.hxxx*(3*xf2_xrd_cu+3*xf_xs2_cu)/6 ... 
                            + derivs.hxxxx*(4*xf3_xs_cu)/24 ... 
                            + derivs.hxxxxx*xf5_cu/120; 
    end 
    % Getting the controls at time t 
    if orderApp==1 
        yVec_cu = model.g0 + derivs.gx*xf_cu; 
    elseif orderApp==2 
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu)+derivs.gxx*(xf2_cu)/2; 
    elseif orderApp==3 
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu+xrd_cu)+derivs.gxx*(xf2_cu+2*xf_xs_cu)/2 ... 
            + derivs.gxxx*(xf3_cu)/6; 
    elseif orderApp==4 
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu+xrd_cu+x4th_cu) ... 
            + derivs.gxx*(xf2_cu+2*xf_xs_cu+2*xf_xrd_cu+xs2_cu)/2 ... 
            + derivs.gxxx*(xf3_cu+3*xf2_xs_cu)/6 ... 
            + derivs.gxxxx*xf4_cu/24; 
    elseif orderApp==5 
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu+xrd_cu+x4th_cu+x5th_cu) ... 
            + derivs.gxx*(xf2_cu+2*xf_xs_cu+2*xf_xrd_cu+2*xf_x4th_cu+xs2_cu+2*xs_xrd_cu)/2 ... 
            + derivs.gxxx*(xf3_cu+3*xf2_xs_cu+3*xf2_xrd_cu+3*xf_xs2_cu)/6 ... 
            + derivs.gxxxx*(xf4_cu+4*xf3_xs_cu)/24 ... 
            + derivs.gxxxxx*xf5_cu/120; 
    end 
    yGrid(k,:) = transpose((yVec_cu-model.g0)); 
    % The controls 
    c_cu = yVec_cu(1,1);
    r_cu = yVec_cu(2,1);
    pai_cu = yVec_cu(3,1);
    w_cu = yVec_cu(4,1);
    l_cu = yVec_cu(5,1);
    output_cu = yVec_cu(6,1);
    mc_cu = yVec_cu(7,1);
    evf_cu = yVec_cu(8,1);
    vf_cu = yVec_cu(9,1);
    p1_cu = yVec_cu(10,1);
    p2_cu = yVec_cu(11,1);
    p3_cu = yVec_cu(12,1);
    p4_cu = yVec_cu(13,1);
    p5_cu = yVec_cu(14,1);
    p6_cu = yVec_cu(15,1);
    p7_cu = yVec_cu(16,1);
    p8_cu = yVec_cu(17,1);
    p9_cu = yVec_cu(18,1);
    p10_cu = yVec_cu(19,1);
    p11_cu = yVec_cu(20,1);
    p12_cu = yVec_cu(21,1);
    p13_cu = yVec_cu(22,1);
    p14_cu = yVec_cu(23,1);
    p15_cu = yVec_cu(24,1);
    p16_cu = yVec_cu(25,1);
    p17_cu = yVec_cu(26,1);
    p18_cu = yVec_cu(27,1);
    p19_cu = yVec_cu(28,1);
    p20_cu = yVec_cu(29,1);
    p21_cu = yVec_cu(30,1);
    p22_cu = yVec_cu(31,1);
    p23_cu = yVec_cu(32,1);
    p24_cu = yVec_cu(33,1);
    p25_cu = yVec_cu(34,1);
    p26_cu = yVec_cu(35,1);
    p27_cu = yVec_cu(36,1);
    p28_cu = yVec_cu(37,1);
    p29_cu = yVec_cu(38,1);
    p30_cu = yVec_cu(39,1);
    p31_cu = yVec_cu(40,1);
    p32_cu = yVec_cu(41,1);
    p33_cu = yVec_cu(42,1);
    p34_cu = yVec_cu(43,1);
    p35_cu = yVec_cu(44,1);
    p36_cu = yVec_cu(45,1);
    p37_cu = yVec_cu(46,1);
    p38_cu = yVec_cu(47,1);
    p39_cu = yVec_cu(48,1);
    p40_cu = yVec_cu(49,1);
 
    f     = zeros(ny+nx,1);
    Euler = zeros(ny+nx,1);
    xf_cup0 = xf_cup;
    for i=1:length(wGrid)
        % Adding shocks to the states 
        xf_cup(1:nx,1) = xf_cup0(1:nx,1) + model.eta*transpose(eGrid(i,:)); 
        % Getting controls at time t+1 
        if orderApp==1 
            yVec_cup  = model.g0 + derivs.gx*xf_cup; 
        elseif orderApp==2 
            xf2_cup   = kron(xf_cup,xf_cup); 
            yVec_cup  = model.g0 + derivs.gx*(xf_cup+xs_cup)+derivs.gxx*(xf2_cup)/2; 
        elseif orderApp==3 
            xf2_cup   = kron(xf_cup,xf_cup); 
            xf3_cup   = kron(xf2_cup,xf_cup); 
            xf_xs_cup = kron(xf_cup,xs_cup); 
            yVec_cup  = model.g0 + derivs.gx*(xf_cup+xs_cup+xrd_cup)+derivs.gxx*(xf2_cup+2*xf_xs_cup)/2 ... 
                + derivs.gxxx*(xf3_cup)/6; 
        elseif orderApp==4 
            xf2_cup    = kron(xf_cup,xf_cup); 
            xf3_cup    = kron(xf2_cup,xf_cup); 
            xf_xs_cup  = kron(xf_cup,xs_cup); 
            xf4_cup    = kron(xf3_cup,xf_cup); 
            xf2_xs_cup = kron(xf_cup,xf_xs_cup); 
            xs2_cup    = kron(xs_cup,xs_cup); 
            xf_xrd_cup = kron(xf_cup,xrd_cup); 
            yVec_cup   = model.g0 + derivs.gx*(xf_cup+xs_cup+xrd_cup+x4th_cup) ... 
                + derivs.gxx*(xf2_cup+2*xf_xs_cup+2*xf_xrd_cup+xs2_cup)/2 ... 
                + derivs.gxxx*(xf3_cup+3*xf2_xs_cup)/6 ... 
                + derivs.gxxxx*xf4_cup/24; 
        elseif orderApp==5 
            xf2_cup    = kron(xf_cup,xf_cup); 
            xf3_cup    = kron(xf2_cup,xf_cup); 
            xf_xs_cup  = kron(xf_cup,xs_cup); 
            xf4_cup    = kron(xf3_cup,xf_cup); 
            xf2_xs_cup = kron(xf_cup,xf_xs_cup); 
            xs2_cup    = kron(xs_cup,xs_cup); 
            xf_xrd_cup = kron(xf_cup,xrd_cup); 
            xf5_cup    = kron(xf4_cup,xf_cup); 
            xf3_xs_cup = kron(xf_cup,xf2_xs_cup); 
            xf_xs2_cup = kron(xf_cup,xs2_cup); 
            xf2_xrd_cup= kron(xf_cup,xf_xrd_cup); 
            xs_xrd_cup = kron(xs_cup,xrd_cup); 
            xf_x4th_cup= kron(xf_cup,x4th_cup); 
            yVec_cup = model.g0 + derivs.gx*(xf_cup+xs_cup+xrd_cup+x4th_cup+x5th_cup) ... 
                + derivs.gxx*(xf2_cup+2*xf_xs_cup+2*xf_xrd_cup+2*xf_x4th_cup+xs2_cup+2*xs_xrd_cup)/2 ... 
                + derivs.gxxx*(xf3_cup+3*xf2_xs_cup+3*xf2_xrd_cup+3*xf_xs2_cup)/6 ... 
                + derivs.gxxxx*(xf4_cup+4*xf3_xs_cup)/24 ... 
                + derivs.gxxxxx*xf5_cup/120; 
        end 
        c_cup = yVec_cup(1,1);
        r_cup = yVec_cup(2,1);
        pai_cup = yVec_cup(3,1);
        w_cup = yVec_cup(4,1);
        l_cup = yVec_cup(5,1);
        output_cup = yVec_cup(6,1);
        mc_cup = yVec_cup(7,1);
        evf_cup = yVec_cup(8,1);
        vf_cup = yVec_cup(9,1);
        p1_cup = yVec_cup(10,1);
        p2_cup = yVec_cup(11,1);
        p3_cup = yVec_cup(12,1);
        p4_cup = yVec_cup(13,1);
        p5_cup = yVec_cup(14,1);
        p6_cup = yVec_cup(15,1);
        p7_cup = yVec_cup(16,1);
        p8_cup = yVec_cup(17,1);
        p9_cup = yVec_cup(18,1);
        p10_cup = yVec_cup(19,1);
        p11_cup = yVec_cup(20,1);
        p12_cup = yVec_cup(21,1);
        p13_cup = yVec_cup(22,1);
        p14_cup = yVec_cup(23,1);
        p15_cup = yVec_cup(24,1);
        p16_cup = yVec_cup(25,1);
        p17_cup = yVec_cup(26,1);
        p18_cup = yVec_cup(27,1);
        p19_cup = yVec_cup(28,1);
        p20_cup = yVec_cup(29,1);
        p21_cup = yVec_cup(30,1);
        p22_cup = yVec_cup(31,1);
        p23_cup = yVec_cup(32,1);
        p24_cup = yVec_cup(33,1);
        p25_cup = yVec_cup(34,1);
        p26_cup = yVec_cup(35,1);
        p27_cup = yVec_cup(36,1);
        p28_cup = yVec_cup(37,1);
        p29_cup = yVec_cup(38,1);
        p30_cup = yVec_cup(39,1);
        p31_cup = yVec_cup(40,1);
        p32_cup = yVec_cup(41,1);
        p33_cup = yVec_cup(42,1);
        p34_cup = yVec_cup(43,1);
        p35_cup = yVec_cup(44,1);
        p36_cup = yVec_cup(45,1);
        p37_cup = yVec_cup(46,1);
        p38_cup = yVec_cup(47,1);
        p39_cup = yVec_cup(48,1);
        p40_cup = yVec_cup(49,1);
 
        % Getting the states at time t+1 
        c_ba1p = model.h0(1,1) + xf_cup(1,1) + xs_cup(1,1) + xrd_cup(1,1) + x4th_cup(1,1) + x5th_cup(1,1);
        muz_cup = model.h0(2,1) + xf_cup(2,1) + xs_cup(2,1) + xrd_cup(2,1) + x4th_cup(2,1) + x5th_cup(2,1);
        d_cup = model.h0(3,1) + xf_cup(3,1) + xs_cup(3,1) + xrd_cup(3,1) + x4th_cup(3,1) + x5th_cup(3,1);
        n_cup = model.h0(4,1) + xf_cup(4,1) + xs_cup(4,1) + xrd_cup(4,1) + x4th_cup(4,1) + x5th_cup(4,1);
        paistar_cup = model.h0(5,1) + xf_cup(5,1) + xs_cup(5,1) + xrd_cup(5,1) + x4th_cup(5,1) + x5th_cup(5,1);
        a_cup = model.h0(6,1) + xf_cup(6,1) + xs_cup(6,1) + xrd_cup(6,1) + x4th_cup(6,1) + x5th_cup(6,1);
        %% Model Specific 
% START DISPLAYING f
f(1,1)=exp(-vf_cu)*(exp(d_cu)*(((exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))/Css^CHI0)^(1 - CHI)/(CHI - 1) - U0d + (PHIzero*exp(n_cu)*(1 - exp(l_cu))^(1 - 1/PHI))/(1/PHI - 1)) - U0 + (AA*BETTA)/exp(evf_cu)^(1/(ALFA - 1))) - 1;
f(2,1)=exp(-w_cu)*(KAPAw*Wss - (Css^CHI0*PHIzero*exp(n_cu)*((exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))/Css^CHI0)^CHI*(KAPAw - 1))/(1 - exp(l_cu))^(1/PHI)) - 1;
f(3,1)=(BETTA*exp(-d_cu)*exp(-pai_cup)*exp(d_cup)*exp(r_cu)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(4,1)=- (exp(-a_cu)*exp(-mc_cu)*exp(w_cu)*exp(l_cu)^THETA)/(Kss^THETA*(THETA - 1)) - 1;
f(5,1)=1 - (exp(-mc_cu)*(ETA + (ZI*exp(pai_cu)*(exp(pai_cu)/PAIss^NU - 1))/PAIss^NU - (BETTA*ZI*exp(-d_cu)*exp(-output_cu)*exp(d_cup)*exp(muz_cup)*exp(output_cup)*exp(pai_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*(exp(pai_cup)/PAIss^NU - 1)*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(PAIss^NU*(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI) - 1))/ETA;
f(6,1)=Rss*exp(PHIl*log(exp(l_cu)/lss) + PHIc*(log(exp(muz_cu)) - log(MUZss) - exp(c_ba1) + exp(c_cu)) + PHIc_1*(log(exp(muz_cup)) - log(MUZss) - exp(c_ba1p) + exp(c_cup)) - PHIpai*(log(exp(paistar_cu)) - log(exp(pai_cu)) + log(PAIss)) - PHIpai_1*(log(exp(paistar_cup)) - log(exp(pai_cup)) + log(PAIss)) + PHIy*log(exp(output_cu)/OUTPUTss) + PHIy_1*log(exp(output_cup)/OUTPUTss)) - exp(r_cu);
f(7,1)=(exp(output_cu)*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1))/(exp(c_cu) + DELTA*Kss) + 1;
f(8,1)=Kss^THETA*exp(-output_cu)*exp(a_cu)*exp(l_cu)^(1 - THETA) - 1;
f(9,1)=exp(-evf_cu)*((exp(vf_cup)*exp(muz_cup)^((CHI - 1)*(CHI0 - 1)))/AA)^(1 - ALFA) - 1;
f(10,1)=exp(-p1_cu)*exp(-r_cu) - 1;
f(11,1)=(BETTA*exp(-d_cu)*exp(-p2_cu)*exp(-pai_cup)*exp(d_cup)*exp(p1_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(12,1)=(BETTA*exp(-d_cu)*exp(-p3_cu)*exp(-pai_cup)*exp(d_cup)*exp(p2_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(13,1)=(BETTA*exp(-d_cu)*exp(-p4_cu)*exp(-pai_cup)*exp(d_cup)*exp(p3_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(14,1)=(BETTA*exp(-d_cu)*exp(-p5_cu)*exp(-pai_cup)*exp(d_cup)*exp(p4_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(15,1)=(BETTA*exp(-d_cu)*exp(-p6_cu)*exp(-pai_cup)*exp(d_cup)*exp(p5_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(16,1)=(BETTA*exp(-d_cu)*exp(-p7_cu)*exp(-pai_cup)*exp(d_cup)*exp(p6_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(17,1)=(BETTA*exp(-d_cu)*exp(-p8_cu)*exp(-pai_cup)*exp(d_cup)*exp(p7_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(18,1)=(BETTA*exp(-d_cu)*exp(-p9_cu)*exp(-pai_cup)*exp(d_cup)*exp(p8_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(19,1)=(BETTA*exp(-d_cu)*exp(-p10_cu)*exp(-pai_cup)*exp(d_cup)*exp(p9_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(20,1)=(BETTA*exp(-d_cu)*exp(-p11_cu)*exp(-pai_cup)*exp(d_cup)*exp(p10_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(21,1)=(BETTA*exp(-d_cu)*exp(-p12_cu)*exp(-pai_cup)*exp(d_cup)*exp(p11_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(22,1)=(BETTA*exp(-d_cu)*exp(-p13_cu)*exp(-pai_cup)*exp(d_cup)*exp(p12_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(23,1)=(BETTA*exp(-d_cu)*exp(-p14_cu)*exp(-pai_cup)*exp(d_cup)*exp(p13_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(24,1)=(BETTA*exp(-d_cu)*exp(-p15_cu)*exp(-pai_cup)*exp(d_cup)*exp(p14_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(25,1)=(BETTA*exp(-d_cu)*exp(-p16_cu)*exp(-pai_cup)*exp(d_cup)*exp(p15_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(26,1)=(BETTA*exp(-d_cu)*exp(-p17_cu)*exp(-pai_cup)*exp(d_cup)*exp(p16_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(27,1)=(BETTA*exp(-d_cu)*exp(-p18_cu)*exp(-pai_cup)*exp(d_cup)*exp(p17_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(28,1)=(BETTA*exp(-d_cu)*exp(-p19_cu)*exp(-pai_cup)*exp(d_cup)*exp(p18_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(29,1)=(BETTA*exp(-d_cu)*exp(-p20_cu)*exp(-pai_cup)*exp(d_cup)*exp(p19_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(30,1)=(BETTA*exp(-d_cu)*exp(-p21_cu)*exp(-pai_cup)*exp(d_cup)*exp(p20_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(31,1)=(BETTA*exp(-d_cu)*exp(-p22_cu)*exp(-pai_cup)*exp(d_cup)*exp(p21_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(32,1)=(BETTA*exp(-d_cu)*exp(-p23_cu)*exp(-pai_cup)*exp(d_cup)*exp(p22_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(33,1)=(BETTA*exp(-d_cu)*exp(-p24_cu)*exp(-pai_cup)*exp(d_cup)*exp(p23_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(34,1)=(BETTA*exp(-d_cu)*exp(-p25_cu)*exp(-pai_cup)*exp(d_cup)*exp(p24_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(35,1)=(BETTA*exp(-d_cu)*exp(-p26_cu)*exp(-pai_cup)*exp(d_cup)*exp(p25_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(36,1)=(BETTA*exp(-d_cu)*exp(-p27_cu)*exp(-pai_cup)*exp(d_cup)*exp(p26_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(37,1)=(BETTA*exp(-d_cu)*exp(-p28_cu)*exp(-pai_cup)*exp(d_cup)*exp(p27_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(38,1)=(BETTA*exp(-d_cu)*exp(-p29_cu)*exp(-pai_cup)*exp(d_cup)*exp(p28_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(39,1)=(BETTA*exp(-d_cu)*exp(-p30_cu)*exp(-pai_cup)*exp(d_cup)*exp(p29_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(40,1)=(BETTA*exp(-d_cu)*exp(-p31_cu)*exp(-pai_cup)*exp(d_cup)*exp(p30_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(41,1)=(BETTA*exp(-d_cu)*exp(-p32_cu)*exp(-pai_cup)*exp(d_cup)*exp(p31_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(42,1)=(BETTA*exp(-d_cu)*exp(-p33_cu)*exp(-pai_cup)*exp(d_cup)*exp(p32_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(43,1)=(BETTA*exp(-d_cu)*exp(-p34_cu)*exp(-pai_cup)*exp(d_cup)*exp(p33_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(44,1)=(BETTA*exp(-d_cu)*exp(-p35_cu)*exp(-pai_cup)*exp(d_cup)*exp(p34_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(45,1)=(BETTA*exp(-d_cu)*exp(-p36_cu)*exp(-pai_cup)*exp(d_cup)*exp(p35_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(46,1)=(BETTA*exp(-d_cu)*exp(-p37_cu)*exp(-pai_cup)*exp(d_cup)*exp(p36_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(47,1)=(BETTA*exp(-d_cu)*exp(-p38_cu)*exp(-pai_cup)*exp(d_cup)*exp(p37_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(48,1)=(BETTA*exp(-d_cu)*exp(-p39_cu)*exp(-pai_cup)*exp(d_cup)*exp(p38_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(49,1)=(BETTA*exp(-d_cu)*exp(-p40_cu)*exp(-pai_cup)*exp(d_cup)*exp(p39_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(50,1)=exp(-c_ba1p)*exp(c_cu) - 1;
f(51,1)=RHOz*log(exp(muz_cu)/MUZss) - log(exp(muz_cup)/MUZss);
f(52,1)=RHOd*log(exp(d_cu)) - log(exp(d_cup));
f(53,1)=RHOn*log(exp(n_cu)) - log(exp(n_cup));
f(54,1)=RHOp*log(exp(paistar_cu)) - log(exp(paistar_cup));
f(55,1)=RHOa*log(exp(a_cu)) - log(exp(a_cup));
% END DISPLAYING f
        Euler(:,1) = Euler(:,1) + wGrid(i,1)*f(:,1); 
    end 
    ResEuler(k,:) = transpose(Euler(:,1));
end 
delete(gcp('nocreate'))
